A procedure for the inversion of f-mode travel times for solar flows 
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ABSTRACT 

We perform a two-dimensional inversion of f-mode travel times to determine near-surface solar flows. 
The inversion is based on optimally localized averaging of travel times. We use finite-wavelength travel- 
time sensitivity functions and a realistic model of the data errors. We find that it is possible to obtain a 
spatial resolution of 2 Mm. The error in the resulting flow estimate ultimately depends on the observation 
time and the number of travel distances used in the inversion. 



Subject headings: Helioseismology, Sun: convection 



Introduction 



In time-distance helioseismology dDuvall et al. 



1993), one measures the wave travel time from one 
point on the surface of the Sun to another. Local flows 
interact with the waves and affect the travel times. 
The travel-time measurements thus contain informa- 
tion about the internal flows. 

One major goal is to infer the near-surface flows at 
the highest spatial and temporal resolution possible, in 
order, for instance, to study supergranulation in detail. 
Previous studies have shown that this is in deed possi 



ble u s ing time-distance helio seismology dGizon et al. 



20001: iDuvall & Gizonll2000l) . To achieve the desired 
accuracy and resolution, a consistent inversion must 
be carried out, which involves two main steps. The 
first step, called the forward problem, is to model the 
wave-flow interaction to obtain sensitivity kernels that 
relate the travel times to small-amplitude flows. The 
second step is to use the kernels and the travel-time 
measurements, together with noise estimates, to infer 



the flows. 

Here we invert f-mode travel times for two - dimen- 
sional, spatially-varyi ng, horizontal flows. We use the 
sensitivity kernels of iJackiewicz et al. (2006), which 
are computed in the first Born approximation. These 
kernels are two dimensional as they are averaged over 
the depth over which the f-mode kinetic energy is sig- 
nificant. We expect the inversion to return horizontal 
flows at an average depth of 1 Mm below the surface. 
We note that the kernels are sensitive to flows at hor- 
izontal scales as small as half the f-mode wavelength 
of about 5 Mm at 3 mHz. 

The inversion scheme adopted in this study is the 
Subtracti ve Optimally Localized Averagin g (SOLA) 
method (Pijpers & Thompson! 1992 . 19941) . which is 
commonly used in helioseismic inversions but has not 
yet been demonstrated for the time-distance technique. 
To keep the presentation simple, we formulate the two- 
dimensional inversion for only one travel distance of 
4.96 Mm. The method can be generalized for many 
distances. 
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In Section [2] we discuss briefly the inputs to 
the inversion, in Section [3] we formulate the two- 
dimensional SOLA inversion scheme, and in Section|4] 
the results of the inversion are shown. We conclude in 
Section [5] 

2. Inputs to the inversion 

We work in a Cartesian coordinate system (x, y) 
representing a small area of the solar surface. The x 
coordinate is measured in the pro-grade direction (to- 
ward the west limb) and the y coordinate is measured 
toward the north pole. The quantities needed for the 
inversion are the travel-time measurements, the sen- 
sitivity kernels, and the noise-covariance matrix. We 
briefly outline these below. 

We choose to consider three different types of f- 
mode travel-time measurements. These are the 'oi', 
'we', and 'ns' travel times, which refer to the out 
minus in, west minus east, and north minus south 
travel-time differences, respectively. These measure- 
ments are very similar to the spatial averages of travel 
times over quadran ts that have been introduced by 
Duvalletal .1 dl997h . A difference is that we are us- 
ing flie_toTCl4ime measurement technique described 
by iGizon & Birch] d2004l) . The 'oi' travel times are 
measured between a center point and a concentric an- 
nulus of radius 4.96 Mm; they are given by the dif- 
ference in travel time between waves that propagate 
outward from and inward toward the center point. The 
'we' travel times are obtained by measuring the travel 
time using the wave signal at the center point and the 
signal averaged over the annulus with a weighting of 
cos 8, where 8 is the angle between the x direction and 
each point on the annulus. Thus, the 'we' travel time 
is mostly sensitive to flows in the x direction. Simi- 
larly, the 'ns' travel time is obtained using a weighting 
of sin 8 to give sensitivity to flows in the y direction. 

We use the short notation T a fa) to denote the 
travel times; the superscript a is one of 'oi', 'we' or 
'ns'. The horizontal vector r*j gives the center position 
of the i-th annulus and is defined on a grid with spacing 
dx — dy = 0.826 Mm in both the x and y directions 
(twice the MDI high-resolution spatial sampling). Let 
us denote the horizontal flow by u = (u x ,u y ). When 
the flow is weak, the relationship between the travel 
times and u can be written as 

= J2 Ka ^ -*■<)• ufa) +* Q (*■<), (i) 
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where the two-dimensional vector-value functions 
K a = (K°,K") are the sensitivity kernels. The 



quantity n a denotes the random error associated with 
the measurement of r a . The sum runs over all spa- 
tial pixels. Note that to simplify the notation, we have 
discretized the flow. 

The f-mode sensitivity kernels K a are a-averages 
o f the point-to-point Bor n sensitivity kernels presented 
in 



Jacki ewicz et al. (2006). The six kernels used in this 



study are shown in Figure [T] 

The statistical properties of the noise in travel time s 
were discussed in detail by IGizon & Birchl (|2004). 
They are specified by the covariance matrix A, 



A^(r % - rj ) = Cov [n Q fa), n^fa)] 



(2) 



We compute A according to equation (28) of Gizon & B irch 
(2004) and using the model f-mode power spectrum 



from Jackiewicz et al 



(2006). x 
elements scale with the observation time T as 1/T. 
To speed up the computation, the full matrix A(r) 
can be constructed using a subset of r values and the 
symmetries of the problem. Figure [2] shows example 
calculations. The diagonal components A aa are gener- 
ally an order of magnitude larger than the off-diagonal 
components A a/3 with a ^ 0. 

All of the inputs to the inversion have now been 
defined and computed. We are now ready to discuss 
the inversion procedure. 

3. SOLA Inversion 

We want to estimate, for example, u x (r ) given a set 
of travel times r Q , the kernels K a , and the noise co- 
variance A af} (Eq. HI). We adopt the SOLA scheme 
developed by Pijpers & Thompson (1992). The strat- 
egy is to search for a set of coefficients w a such that 
an estimate of u x {r) is 

t£(r) := wa ( r i ~ r ) T "fa) « u *(r). (3) 



The above definition of u x is equivalent to 

Thi(r) = ^K{r k ~r)-u{r k )+^w a {r j -r)n a {r J ), 

k j,a 

(4) 



where 



K(r):=^ a (rj)K a (r- rj ) (5) 



is the averaging kernel. The averaging kernel is vector 
valued, i.e. 7C = (JC Xl K. y ). Ideally the x-component 
of the averaging kernel, K, x , should be a delta function 
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Fig. 1. — Travel-time sensitivity kernels used in this study, K a = (K x , K£), where a stands for 'oi\ 'we', 'ns'. The 
units of the kernels are s (m/s) -1 . The grid spacing is dx = dy = 0.826 Mm. The center of each annulus lies at 
r = (0, 0). The spatial sum of K x c and Ky S is —0.043 s (m/s) -1 . The sums of all other kernels are zero. 



of position and JC y should vanish in order to recover 
u x perfectly. This is not possible in general because 
of noise and a limited set of travel times. Instead, in 
the spirit of SOLA, we attempt to choose weights so 
that JC x (r) will resemble a prescribed target function, 
T x (r), which is localized around r = 0. In practice, 
we take a Gaussian function for T x (r ) with dispersion 
cr, such that 



T(r) 



-r 2 /2a 2 



(6) 



where r 



2vrcr 2 

and || • || is the 2D vector norm. The 
total integral of T x (r) is one. If we can construct an 
averaging kernel K w T with a small, then u^(r) will 
be an estimate of u x {r') in the neighborhood \\r' — 
r|| < a. As a result we call cr, or FWHM= 2er\/21n 2, 
the target resolution. The effective resolution depends 
on how well the averaging kernel matches the target 
function. 

In general, OLA inversions seek a balance between 
the resolution of the solution and the error magnifica- 



tion that propagates through the inversion due to the 
uncertainties of the measurement. To quantify this 
more precisely, we define two measures. The first mea- 
sure gives the misfit between the averaging kernel and 
the chosen target function, 



misfit 



\K(n) - T(n) 



(7) 



One wants the misfit to be as small as possible. The 
second measure quantifies the error in u^: 



w a {r i )K a l i {r i -r j )w' i {r j ). (8) 



It is also desirable to keep this term as small as possi- 
ble. The balance of the two measures defined in equa- 
tions © and dS) is achieved by minimizing 



misfit + /i error 2 



(9) 



with respect to w a , where /j, is a trade-off, or regular- 
ization parameter. The minimization of expression (O 
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Fig. 2. — Example calculations of the noise covari- 
ance, A Q/3 . The left column shows the maps A OIOI (r) 
(top), A oiwc (r) (middle), and A wowc (r) (bottom) as 
functions of r = (x,y). The units are s 2 and the ob- 
servation time is T = 12 hr. The right panels are cuts 
through the corresponding maps at y = 0. 



is also subject to the constraint 



5>(r, 



(10) 



The trade-off parameter is varied through a range of 
values. The value for /i is chosen so that one is content 
with the amount of misfit and the amount of error that 
results from the inversion. 

Upon carrying out the minimization of expres- 
sion (0 with respect to the weights, one finds that 
for all indices k and measurement types 7, 



E 



A k ~ / j ,w a (r j ) = fcfc 7 , 



(11) 



where we have defined 

A kl , ja := ^^(ri-r.O-^Cn-rfe) 

i 

+ LtA~< a (r k -r j ), (12) 
b kl := ^^(ri-rfcJ-TCr,). (13) 

i 

In addition, the constraint (fTTTb implies 



E 



E*>*-^') 



w a (r 3 ) = x 



(14) 



Taken together, equations (fTTT i and (TBI lead to a sys- 
tem of linear equations that can be written as 



Aw; = b. 



(15) 



A computational advantage of the SOLA formulation 
is that the target function is contained in b, and so 
the matrix A need not be inverte d for each new target 
that s atisfies the same constraint (IPiipers & Thompson 
19921) . 

There are inherent problems in the solution given 
by A _1 b, since the matrix A may be ill-conditioned. 
The matrix A may have small singular values which 
result in oscillatory singular vectors: upon inversion 
of the matrix, the high-frequency components are am- 
plified and the corresponding solution tends to be 
meaningless. This is a commonly occurring prob- 
lem in inversions. To rectify this, we regularize the 
two matrices corresponding to the two terms in equa- 
tion ( fT2l by performing a singular value decomposi- 
tion (SVD) and removing the smallest singular values, 
a process referred to as truncated SVD. This tech- 
nique is commonly used in helioseism ic inversions 
(e.g JChristensen-Dalsgaard et al ] |l993h . 



4. Results 

We solve equation (TT~5T > for many sets of inversion 
coefficients w a using a range of values for the trade-off 
parameter /1. This is done for different FWHM values 
of the target function as well. The results can be stud- 
ied in error-m isfit space by computing the 'L-curves' 
( Hansenfl998h . shown in Figure [3] 

Each curve in Figure [3] corresponds to a particular 
value of the FWHM, and each point on each curve cor- 
responds to a different set of weights. In this case, 
there are 45 values of /1 which span over 15 orders of 
magnitude to produce each L-curve. The goal is to 
choose a particular set of weights so that one is satis- 
fied with the trade-off between the misfit and error. 



4 



10' 



10 



10 



10' 



1^=2.156-09 s m 
(i = le-05 s m 



4.96 Mm 
5.78 Mm 




10 



10 ' 

misfit 2 (arb. units) 



10 



P 



0.02 

0.015 -g 
0.01 t. 

: 0.005 





Avg. Kernel 




i 









11.02 
0.015 
0.01 

0.005 




error= 38.5 m/s / 


\ Avg. Kernel 




error= 25.6 ni/s 


, , Avg. Kernel 


misfit= 7.1e-05 / 




0.02 


misfit= 0.039 


,' \ Target 






0.015 






"i / 




0.01 


^2 








0.005 







Fig. 3. — A set of L-curves from the inversion. Each 
curve represents a different FWHM of the target func- 
tion, noted at the top left of the curve, and each point 
corresponds to a different value of the regularization 
parameter /i. The circles and squares denote where /Ui 
and H2 lie on each curve, respectively. 

Highlighted in Figure [3] are two specific values 
of the regularization parameter, /j, = fii and /i = 
/i2, which represent two possible trade-off choices. 
The points on the L-curves given by H> = H2 are 
ones chosen by an algorithm that finds the 'corner' 
of discreet L-curve s for general inverse problems (see 
Hansen et al. 2006). The 'corner' is generally consid- 
ered to be the point of the best trade-off. However, in 
this example, choosing the weights at the points given 
by n = fXi seems much more appropriate. For in- 
stance, consider the L-curves with the larger values of 
FWHM in Figure [3] As one increases fi from /ii to 
the increase in misfit is nearly six orders of magni- 
tude, while the decrease in error is only about a factor 
of two. 

The difference between hi and fi2 is further ex- 
plored in Figure |U where we examine the misfit more 
closely by plotting cuts through the averaging ker- 
nels JC X and the target function T x , at hi and /j,2 and 
FWHM=4.96 Mm. For this example JC y — accord- 
ing to equation fllOl . One sees that even though the 
error at /i = \ii is modestly lower than the error at 
/i = hi, the averaging kernels at fi = fj,2 do not 
match the target function very well. In particular, there 
is some structure away from the central peak, which 
is undesirable, while the averaging kernel and target 



Fig. 4. — Averaging kernels JC X for two different val- 
ues of h, at fixed FWHM= 4.96 Mm. The 2 panels 
on the left are for (ii> while the two on the right are 
for fi2- The top row is the averaging kernel for each 
fjL. The bottom row shows a cut along the y — line 
for each averaging kernel (solid line), as well as a cut 
through the target function (dashed line). The value of 
the error and misfit is also given for each /i. The spatial 



sum of each JC T is 1. T=12 hrs. Note that JC h 



0. 



match almost exactly at fi = Therefore, in what 
follows, we choose to study the set of weights corre- 
sponding to h = /ii for each FWHM. With this choice 
of vanishing misfit, one may use the terms resolution 
and FWHM interchangeably. 

We note, however, that for a very small FWHM, as 
in Figure|5] it might be reasonable to choose a slightly 
larger value for /.i than fj, = fj.1- The L-curve for this 
FWHM =1.65 Mm (Fig. [3]) demonstrates that the error 
decreases very rapidly with very little increase in misfit 
as h is increased. However, because we are using only 
one travel distance in this study, we will not consider 
such small FWHM any further since the error is too 
large. Using many travel distances would be necessary 
to lower the noise at such a high resolution. 

Another helpful way to visualize the trade-off 
is to study plots of contours of constant error and 
misfit in the FWHM - u space, as in Figure [6] 
(Piipe rs & Thompson! 1994 ). This figure shows again 
that if one demands a small misfit, as in the middle to 
upper left region of the plot, achieving a small error 
requires settling for a poorer resolution. Ideally, one 
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Fig. 5.— Same as Figure g] but for FWHM 
1.65 Mm. 




Fig. 6. — Contour plot of error and misfit from 
an inversion, as functions of FWHM and log 10 fi. 
The solid lines are contours of constant error, in 
units of m s _1 . The dotted lines are contours of 
constant log 10 (misfit), whose values are given in 
the boxes. The left and right stars correspond to 
FWHM=4.96 Mm and [ii and H2> respectively. As in 
all previous plots, T = 12 hrs. 



would like to be in the lower right part of the plot, 
where a small FWHM (high resolution) accompanies 
a small error. However, the misfit there is large and the 
averaging kernels are not localized at all. 

The inversion coefficients w a corresponding to fi = 
Hi and FWHM = 4.96 Mm are given in Figure|7] Ac- 
cording to equation J3J, the weights average the travel- 
time measurements to give an estimate of the solution. 

The typical dependence of the error on the FWHM 
for an observation time T=12 hrs is shown in Figure[8] 
In order to achieve a reasonable error (say < 20 m s _1 ) 
for T=12 hrs, one cannot hope for sub-wavelength res- 
olution (< 5 Mm). Therefore, to increase the reso- 
lution and keep the same error estimate, it becomes 
necessary to observe for a longer time. As has al- 
ready been mentioned, the covariance of the random 
error in the travel-time measurements scales as T _1 
(see Gizon & Birch 2004). We use this result to esti- 
mate the errors over a wide range of typical observa- 
tion times, as shown in Figure|9] The figure shows that 
to reach sub-wavelength resolution with a reasonable 
error, the observation time should span several days. 
Conversely, Figure [9] is also very useful for choosing 
values of T and resolution that will give an acceptable 
noise level. 

5. Discussion 

We have formulated and solved a two-dimensional 
SOLA inversion of f-mode travel times for flows. It 
was shown that this inversion procedure works, in 
that there exist averaging kernels that match the tar- 
get function very well. The random noise on the flow 
estimates, however, has been found to be quite large 
for one travel distance and typical values of the obser- 
vation time, T. To lower the noise, it will be necessary 
to use many travel distances. 

We suggest that the inversion coefficients should be 
chosen to give essentially the smallest misfit between 
the target function and the averaging kernel. A use- 
ful consequence of this particular choice is that one 
immediately knows the resolution. Another advantage 
is that travel-time measurements can be inverted sep- 
arately for each annulus radius at fixed target resolu- 
tion. This would allow for a straightforward averaging 
of the results of the inversions for the various annulus 
radii. 

How should one know which target resolution to 
choose? A possible strategy is to choose first an ac- 
ceptable level of random error on the flow estimates 
and to then choose the target resolution and observing 
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Fig. 8. — Error as a function of the FWHM of the tar- 
get function. The two dashed vertical lines denote the 
value of one wavelength and one-half wavelength of f 
modes at 3 mHz. 
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Fig. 9. — Contour plot of error for different observing 
times T and FWHM. The lines of constant error are in 
units of m s _1 . The error is computed using the set 
of inversion coefficients at [i = fii. The two dashed 
vertical lines denote the value of one wavelength and 
one-half wavelength of f modes at 3 mHz. 



time accordingly (see Fig.|9]l. 

We note that the present inversion procedure can 
presumably easily be generalized to infer flows in three 
dimensions (the size of the matrix to be inverted would 
remain the same). The generalization to other types of 
solar perturbations besides flows is straightforward, as 
long as the corresponding travel-time sensitivity ker- 
nels are available. 
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Fig. 7. — Inversion coefficients w a , corresponding to fi = FWHM=4.96 Mm, and T = 12 hrs. The units are 
(m/s) s _1 . The averaging kernel is by definition obtained by convolving these coefficients with the kernels according 
to equation (0. 
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